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Abstract 


The  application  of  conformal  mapping  methods  to 
the  solution  of  free-surface  flow  problems  is  considered. 
Methods  of  numerical  conformal  mapping  based  on  Fourier 
series  are  extended  to  handle  efficiently  problems  with 
time-dependent  boundaries.  They  are  shown  to  be 
practicable  only  for  moderately  distorted  geometries. 
Extensions  of  the  Menikof f-Zemach  method  to  'breaking' 
geometries  are  presented.  These  latter  methods  are 
robust  at  quite  large  distortions,  but  degrade 
prematurely  in  time-dependent  problems  at  amplitudes 
smaller  than  achieved  by  our  recent  vortex  methods. 
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INTRODUCTION 


In  this  paper,  we  investigate  the  application  of 

conformal  mapping  to  the  solution  of  time -dependent 

potential  flow  problems,  such  as  Rayleigh-Taylor  instability 

and  water  waves.  We  begin  by  formulating  the  water  wave 

problem.  For  two-dimensional,  incompressible,  irrotational , 

free-surface  flow,  the  vclcr  ity  is  expressible  as  v  =  V <J> , 

where  the  potential  4>  satisfies  Laplace's  equation 
2 

V  =  0  in  the  region  y<  n(x,t)  beneath  the  free 
surface  y  =  n(x,t).  Since  the  free  surface  moves  with 
the  fluid, 

=  V  it  =  *y  /  (y  =  n(x't)}'  (1*1) 

where  D/Dt  is  a  Lagrangian  derivative.  Bernoulli's 
law  is  satisfied  throughout  the  fluid  so  that 

§|-  =  -  n(x,t)  +  y  -  pg  (1.2) 

at  the  free  surface  y  =  n(x,t),  where  the  gravitational 
acceleration  is  normalized  to  unity  and  pg  is  the 
applied  surface  pressure.  It  is  assumed  below  that  the 
free  surface  is  periodic  in  x  with  wavelength  2ir. 
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In  order  to  march  forward  in  time,  it  is  necessary 
to  know  V<p  at  the  free  surface.  If  $  is  known  then 
its  tangential  derivative  3<{>/3s  is  computable  but  its 
normal  derivative  3<J)/3n  must  be  found  by  solving 
Laplace's  equation.  Green's  third  formula  expresses 
3<J>/3n  in  terms  of  <p  : 


I  ~~  £n|p-q|<Mq)dq  -  / 
3D  dn  ~  ~  ~  ~  3D 


Sin  |p-q|  |^(q)dq  =  2fT<f>  (p)  . 

(1. 


Here  p,q  are  vectors  lying  on  the  boundary  3D  of  the 

region  D.  Eq.  (1.3)  is  a  linear  integral  equation  of 

the  first  kind  for  the  unknown  function  3<j>/3n.  Once 

3<J>  /3n  has  been  calculated,  Eq.  (1.1)  and  (1.2)  may 

be  used  to  update  the  free  surface  and  potential. 

Numerical  solution  of  (1.3)  for  3<J>/3n  involves 

the  approximation  of  its  logarithmic  kernel  by  a  finite 

matrix.  If  the  continuous  boundary  is  approximated  by  N 

discrete  points,  the  ODeration  count  for  the  solution  of 

3 

the  resulting  linear  system  is  0 (N  )  since  the  matrix 

is  full.  In  addition,  storage  of  the  matrix  requires 
2 

0 (N  )  memory  locations v  For  large  N  the  computational 


costs  are  prohibitive. 


Our  work  is  motivated  by  the  desire  to  develop 

algorithms  with  decreased  operation  counts  and  storage 

requirements  for  solving  free  surface  potential  flow 

problems.  We  have  recently  proposed  a  new  vortex 

1,2 

method  to  solve  these  problems  with  0(N)  memory  and 
2 

0 (N  )  operations  per  time  step.  In  the  present  paper, 

two  conformal  mapping  methods  are  studied.  Both  require 

only  0 (N)  memory.  The  methods  discussed  in  Sec.  2 

require  only  0  (N  X-og ^  N)  operations  per  time  step  but 

are  effectively  limited  to  modest  surface  deformations. 

In  Sec.  3,  modifications  of  the  Menikof f-Zemach  method^ 

2 

that  require  0 (N  )  operations  per  time  step  are  introduced. 
Larger  surface  deformations  can  be  handled  accurately  by 
these  latter  methods. 

2.  CONFORMAL  MAPPING  USING  FOURIER  SERIES 

In  this  Section,  numerical  methods  are  developed 
to  compute  the  conformal  map  z(c)  of  the  unit  disk 
(c|<l  onto  a  simply  connected  finite  region  D  in  the 
complex-z  plane.  A  map  of  the  unit  disk  onto  a  semi¬ 
infinite  periodic  region  R:  w  =  x+iy,  0  £  x  2tt,  y<  n  (x)  , 
is  given  by 

w  =  i  £n  z  (C)  (2.1) 

where  z(£)  is  a  map  of  the  unit  disk  onto  the  interior 
of  the  region  with  boundary 


z  =  exp  t-ix  +  n (x) ] . 


This  sequence  of  conformal  maps  is  depicted  in  Fig.  1. 

Before  proceeding  to  the  discussion  of  methods  to 
compute  the  conformal  map  z{£),  we  note  that  knowledge 
of  z(0  allows  efficient  solution  of  potential 
problems  in  the  region  R.  If  <J> (w)  is  harmonic  in  R 
then  <J>(i  iln  zU))  is  harmonic  in  the  unit  disk. 
Therefore,  the  Dirichlet  problem  can  be  solved  by 
Poisson's  formula.  Also,  since  conformal  maps  are  angle 
preserving,  the  normal  derivative  3<j>/3n  of  cp  on  the 
boundary  ;  R  is  related  to  the  radial  derivative 
3<J>/3p  of  <f>  on  |?|  =  p  =1 : 


y=n (x) 


=  H  2(0 

3p  P=1  dz 


(2.2) 


d?  P  =1 


The  derivative  dz/d?  can  not  vanish  for  |  £  |  <L  if 
z(S)  is  single-valued. 

Let  us  begin  by  characterizing  the  analytic  character 
of  z(C)  in  terms  of  Fourier  series.  The  boundary  values 
z(em)  of  the  conformal  map  z  ( e; )  are  a  periodic  function 
of  the  angle  a  on  the  unit  disk  so 


T 
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z(e10)  =  l  A  elkCl  (2.3) 

k=-®  K 

The  condition  that  z(£)  be  analytic  is 

Ak  =  0,  k  <  0  (2.4) 

and,  in  this  case,  z(s)  is  given  explicitly  by 
00  v 

z(U  =  l  A  (2.5) 

k=0  K 


In  other  words,  an  analytic  transformation  of  the  unit 
disk  onto  a  region  D  is  equivalent  to  a  parametrization 
of  3D  in  terms  of  a  such  that  the  Fourier  representation 
of  3D  has  only  positive  frequency  components. 

Now  we  consider  a  discrete  approximation  to  the 
conformal  map.  Consider  the  equally  spaced  discrete  points 

oij  =  oj  (j  =  0,  ...  N-l)  ,  where  o  =  2tt/n,  and  the  associated 
points  Zj  on  3D.  Then  z^  can  be  represented  as  the 
finite  Fourier  series 


z.  =  z(eia^)  =  l  a.  eika^  (0£j<N) . 
J  N  N  * 


(2.6) 


It  may  easily  be  shown  that 
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t 


a 


k 


l 

P?0 


Ak+pN 


(2.7) 


9.  m  *%.  m  m 

wucr  wa^ 


to  determine  an  approximation  to  the  conformal 


map  z(£)  is  to  require  that  a^  =  0  for  k<  0.  Indeed, 
if  N  is  large  enough  that  is  negligible  for  k>N, 

then  a^  is  negligibly  small  for  k<  0.  This  idea  may 
be  used  to  obtain  iterative  methods^-**  based  on  the  fast 
Fourier  transform  (FFT)  to  compute  the  approximate  conformal 
map.  These  methods  typically  require  0 (N  £og2N)  operations 
per  iteration.  We  note  that  for  any  satisfying 

a^  =  0  for  k  <  0,  the  resulting  conformal  map 


z(C)  =  l  a  £k 
0<k<N/2  * 


(2.8) 


satisfies  z(e10^)  =  z.. 

1 

N  equally  spaced  points 
3D. 


Thus , 
eia^ 


the  map  z  transforms  the 
into  points  z^  lying  on 


As  an  alternative  to  these  iterative  methods,  we  have 
obtained  a  differential  equation  which  relates  the  time  rate 
of  change  of  the  conformal  map  to  the  time  rate  of  change 
of  a  moving  boundary.  This  differential  equation  is 
well  suited  to  the  solution  of  free  surface  flow  problems 
where  the  solution  of  the  potential  problem  determines  the 
time  rate  of  change  of  the  free  surface. 
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Let  the  boundary  be  represented  for  all  time  t 
by  the  equation 


F  (z  (a ,  t)  ,  z(a,t),t)  =  0.  (2.9) 


Differentiation  of  (2.8)  with  respect  to  t  yields 


dF_  3F^  9 z  8F  8z 

at  az  at  az  at 


(2.10) 


and  differentiation  with  respect  to  the  angle  a  (see  Fig.  1) 
gives 


3F  3z  (  3F  9z 
az  3a  3^  8a 


(2.11) 


The  relation 


3F 

1  at 

2  9F  9z 
a  z  aa 


(2.12) 


is  obtained  by  substituting  (2.11)  into  (2.10).  This 

9  z  9  z 

equation  only  provides  the  imaginary  part  of  -^rZ-K— 

a  t  d  CX 

3  Z  3  Z 

The  real  part  of  r-p/^-  is  determined  by  requiring  that  it 
be  analytic  in  the  domain  described  by  (2.9). 

The  right  hand  side  of  (2.12)  is  real  and  can  be 
represented  by  a  conjugate  symmetric  Fourier  series: 
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i  !E 

I  3t 


3F  3_z 
c«z  da 


=  Re 


|n-1 

I  b  exp(ima) 
m=0  m 


Therefore,  analytic  continuation  of  (2.12)  gives 


3  Z  ^  3  Z  r 

3t  “  3a  ^  b  exp(ima) 
m=0  m 


(2.13) 


where 


3z 

3a 


I 

k=0 


(ik)ak  exp  ika) 


(2.14) 


The  right  hand  side  of  Eq.  (2.12)  is  directly  related 
to  the  normal  velocity  of  the  moving  boundary.  In  Cartesian 
coordinates 


F(z,z,t)  =  y  -  n(x,t) 


(2.15) 


so 


3F  =  3n 
3t  ~  3 1 


(2.16) 


Also, 


3F  3_z_ 
9z  3a 


1  ,3F 

2  l3x 


i 


3F.  , 3x 
3y'  l3a 


+ 


i 


(2.17) 
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Substitution  of  (2.16)  and  (2.17)  into  (2.12)  gives 


Im 


in 

at 

3x  3jt_  a^ 
3a  3x  aa 


(2.18) 


The  free  surface  condition  (1.1)  can  be  written 


as 


a n  _  a <()  a n  34> 
at  ~  3y  3x  3x 


(2.19) 


Using  the  decomposition  of  a<J>/3x  and  a4>/3y  into 
tangential  and  normal  components 


a<f> 

fix 

a<t> 

3y 

\  as 

as 

a<t> 

/  ax 

3<t  . 

37  - 

\  as 

as 

(2.18) 

gives  the 

Ts  3n 


) 


(2.20) 


coordinates 


Im 


3z 

at 

az 

3a 


p  and  a  s 


_  3  <p 
~  ~  Bp 


z 

3z 

p=i 

3a 

(2.21) 


An  analogous  result  was  obtained  using  perturbation  methods 

O 

by  Kantorovich  and  Krylov. 
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Eq.  (2.21)  describes  the  motion  of  points  following 
the  conformal  map  of  the  free  surface  rather  than  the 
Lagrangian  or  Eulerian  paths.  Bernoulli's  equation 
(1.2)  must  also  be  modified  to  take  this  fact  into 
account : 


51  = 

Dt 

!  -  Jin 

!  zl 

1  1  /H\  2  -34) .  2 

2  v3pJ  l3a; 

z 

3  z 

L_  _ 

3a 

+  il 

Dx 

+  il  5Z 

3x 

Dt 

3y  Dt 

Once  the  conformal  map  (2.4)  from  the  unit  circle 
is  known,  the  solution  of  the  Dirichlet  problem  may  be 
given  in  terms  of  a  Fourier  series 


<J>  (pela) 


l 

n=0 


[g 


n 


cos  na  +  h  pn  sin  n  a] 


n 


n 


(2.23) 


On  the  unit  circle  (p  =  1)  the  tangential  and  normal 
derivatives  are  given  by 


3  <J> 
3a 


34> 

9p 


n=0 


n=0 


[-ng  sin  na  +  n  h  cos  na  ] 
3n  n 


[ng  cos n a  +  nh  sin n  a  ]  , 
^n  n 


(2.24) 


(2.25) 


which  are  computable  using  FFT's  in  0 (N  £og2N)  operations 
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Now  let  us  summarize  the  steps  involved  in  marching 
from  time  t  to  t  +  At  by  this  method.  At  time  t, 
the  points  {z^}  and  potentials  { 4> ^ }  are  assumed  known. 
First,  the  coefficients  gn,  hR  in  (2.23)  are  obtained 
from  using  an  FFT .  Next,  ^(j>  is  computed  on  the 

boundary  using  (2.24)  and  (2.25).  Then  <f> is 
updated  by  (2.22)  and  z^  is  updated  by  (2.21). 

The  total  operation  count  is  0(N  £o^N)  .  Note  that  the 
conformal  map  is  uniquely  defined  by  (2.21)  and  the 
supplementary  conditions  a^  =  0,  Im(a^)  =  0. 

As  a  test  of  this  time  dependent  mapping  method, 
we  study  the  propagation  of  Stokes'  permanent  water  waves. 
In  a  frame  of  reference  moving  with  the  wave  speed,  the 
numerically  calculated  profiles  should  be  steady.  The 

initial  conditions  for  the  calculation  are  obtained  using 
Pade  approximants  of  perturbation  expansions  of  the  Stokes 
waves11.  The  time-dependent  equations  (2.21),  (2.22) 
are  solved  by  a  fourth-order  Adams-Moulton  predictor 

corrector  scheme.  As  with  other  simulations  of  propagating 

2  9 

nonlinear  water  waves,  '  an  instability  of  the  free 

surface  quickly  develops  unless  damping  is  applied.  To 
remove  this  instability,  we  periodically  apply  a  five- 

.  9 

point  smoothing  operator  . 

For  a  Stokes  wave  with  peak-to-trough  amplitude  80% 
of  the  maximum  allowed  by  theory,  we  choose  the  time  step 
to  be  2tt/400  with  wavelength  2tt  and  apply  smoothing 
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every  tenth  step.  The  resulting  wave  profile  is  plotted 
in  Fig.  2.  The  dots  indicate  the  position  of  the  points 
Zj  used  to  calculate  the  conformal  map  of  the  wave. 

The  solid  line  is  the  wave  profile  computed  by  Pade  approximants 
translated  by  an  amount  equal  to  the  nonlinear  phase  speed 
multiplied  by  the  time.  The  computation  time  for  one 
evaluation  of  the  time  derivatives  of  the  map  and  potential 
is  3ms  on  the  CRAY-1  computer  using  N  =  64  points.  The 
total  computation  time  for  the  motion  of  the  wave  through 
one  period  is  about  Is. 

The  conformal  mapping  method  described  here  works 

well  provided  the  region  is  not  highly  distorted.  As 

the  region  of  interest  becomes  more  distorted  the  points 

corresponding  to  the  conformal  map  tend  to  crowd3,10. 

For  example,  consider  the  conformal  mapping  from  the  unit 

disk  to  the  region  lying  below  y  =  A  coskx.  The  number 

of  terms  N  that  must  be  retained  in  the  Fourier  expansion 

(2.4)  to  obtain  a  good  representation  of  this  map  satisfies 

Jin  N  ~ #kA  as  kA  ->■«>, 3  where  X  =  \  f  sl^-X  dx  -  2,909, 

z  0  x 

showing  the  difficulty  of  mapping  from  the  unit  disk  to  a 
deformed  region. 

When  N  is  large,  almost  all  of  the  equally  spaced 
points  e  a 3  on  the  unit  circle  are  mapped  into  points 
Zj  that  are  crowded  into  small  intervals  on, the  boundary 
of  the  domain  D.  Dubiner  has  recently  made  a  detailed 
analysis  of  this  problem  and  has  shown  that  the  crowding 
occurs  whenever  the  region  being  mapped  has  a  'narrow' 


section.  This  effect  occurs  in  high  amplitude  Rayleigh- 
Taylor  instability  and  in  breaking  waves.  The  FFT  method 
is  not  effective  in  dealing  with  these  highly  distorted 
geometries . 

When  the  domain  D  is  highly  deformed,  the  iterative 
4-8 

methods  and  our  differential  equation  method  do  give  a 
conformal  map  of  the  unit  disk  onto  a  domain  that  passes 
through  the  desired  points  of  3D.  However,  unless  N 

is  unreasonably  large,  the  conformal  map  so  obtained 
will  have  large  deviations  from  3D  between  the  points 
z j .  Indeed,  (2.8)  gives  an  accurate  conformal  map  of 
3D  only  if  a^  decreases  rapidly  as  k  increases  to 


One  possible  approach  to  the  crowding  problem  is 

to  use  a  sequence  of  mappings  of  the  disk  onto  successively 

more  highly  deformed  regions.  Such  iterated  mappings 

are  still  under  study.  For  such  methods,  one  result  seems 

assured,  namely  that  the  operation  counts  must  degrade 

2 

from  0(N  ^0^2  N)  to  0 (N  )  or  worse.  In  this  case, 
these  Fourier  series  methods  are  probably  inferior  to  the 
methods  to  be  described  in  Sec.  3. 

3.  APPLICATION  OF  THE  MENIKOFF-ZEMACH  METHOD 

The  Fourier  series  methods  for  mapping  the  unit  disk 
onto  D  can  not  accurately  handle  highly  distorted  domains 
as  the  crowding  phenomenon  causes  a  severe  loss  of  resolution 
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m  some  part  of  the  physical  boundary.  This  difficulty 
may  be  overcome  by  mapping  D  onto  the  unit  disk  with 
a  regular  distribution  of  points  on  3d.  The  crowding 
then  occurs  on  the  boundary  of  the  unit  disk.  Even 
with  a  highly  nonuniform  distribution  of  points  on  the 
unit  circle,  the  potential  problem  in  the  unit  disk  is 
still  readily  solved  by  Poisson's  formula. 

Recently,  Menikoff  and  Zemach^  have  developed  a 
new  nonlinear  integral  equation  for  conformal  mapping 
of  the  region  R  above  y  »  n (x)  onto  the  periodic 
semi-infinite  strip  S:0<u<2tt,  0<v<«>  .  Their  method 

requires  relatively  few  points  to  achieve  accurate 
results  for  distorted  domains. 

A  simple  extension  of  Menikoff  and  Zemach's  equation 
which  is  valid  for  general  periodic  interfaces  is 
derived  here  and  is  used  to  investigate  the  crowding 
phenomenon  for  multivalued  (or  'breaking  wave')  interfaces. 
A  time  dependent  version  of  the  equation  is  also  developed. 
This  approach  reduces  to  the  integration  of  N  nonlinear 
differential  equations. 

The  Menikof f-Zemach  equations,  generalized  to  handle 
conformal  maps  of  a  domain  with  boundary  curve  parametrized 
as  x  =  x  (e)  ,  y  =  y  (e)  ,  are 


3 


y  (e) 


+  2 


2ir 

/  in 
0 


sin  2 (u (e) -u  (e  * ) ) 


sm 


J(e-e') 


dx 

de' 


<e,i  %r 


2tt 

+  /  cot 

0 


( 


e-e 

2 


I 


[(x(e')-e')  -  (x(e)-e)] 


(3.1a) 


u  (e)  =  x  (e) 


x  +  2 

oo 


2tt 


/  in 

0 


1 

sin  2 (u (e) -u (e ' ) ) 
sin  j(e-e ' ) 


dy  de ' 
de '  2ir 


+ 


2ir  , 

/  cot  (S=e_) 
0  ^ 


(y (e' ) -y (e) ] 


de' 

2n  * 


(3.1b) 


Here  e  is  chosen  so  that  x(0)  =  0,  x(2ir)  =  2tt,  u(e) 
is  defined  so  (x(e),y(e))  is  mapped  into  (u(e),0),  and 
y<x>'xco  are  determined  by  the  condition  that  u(0)  =  0. 

Note  that  (3.1a)  and  (3.1b)  are  equivalent;  either  one 
can  be  used  to  determine  u(e).  Once  u(e)  is  found  by 
solution  of  (3.1),  the  conformal  map  is  determined. 

Eqs.  (3.1)  are  derived  from  the  pair  of  Hilbert 
transforms : 

2tt  ,  .  , 

Re  t  G  ( u ,  0 )  ]  =  Re(Gj  -  P  /  cot  (~^~)  Im  [G  (u '  ,  0]  ~ 


(3.2a) 
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Im[G  (u,  0) 


Im(G  )  +  P 

CO 


r  cot  )  Re  [G  (u '  ,0)  ] 


(3.2b) 


where  G(W)  is  analytic  in  the  upper-half  W-plane  and 
G (w)  =  G^  +  0(l/|w!)  as  jw  j  -*■“>.  Eq.  (3.1)  follows 
if  G  =  x  +  iy-w  ,  where  W  =  u+iv. 

The  Hilbert  transforms  (3.2)  are  also  useful  for 
solving  potential  problems  in  the  region  R.  If  the  map 
function  u(e)  is  known,  boundary  values  of  a  potential 
4>  on  3R  may  be  related  to  corresponding  boundary  values 
of  a  potential  $  defined  in  the  strip  S  in  the  W-piane: 


4>  (u)  =4>  (e) 


34>  /du 
3s/ de 


3jM  _  3j£  ds  /du 
3v.L_q  3n  de  /  de 


(3.3) 


where  s  and  n  are  the  tangential  and  normal  directions 
to  3R.  The  tangential  and  normal  derivatives  of  $  are 
the  real  and  imaginary  parts  of  an  analytic  function  in  S 
so  they  are  related  by  the  Hilbert  transform: 


3<t>  ds 
3n  de 


2tt 

/  cot 
0 


(3.4) 


•» 


-17- 


r 
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Note  that  in  the  application  of  (3.4)  it  is  necessary 
to  compute  du/de  with  some  care.  We  have  found  it  best 
to  find  du/de  by  using  the  Hilbert  transform  of 
Jin  dz/dw  to  obtain  an  equation  for  Jin  du/de. 

In  order  to  examine  the  crowding  properties  of  domains 
bounded  by  breaking  waves,  we  use  (3.1)  to  compute  the 
function  u(e)  for  the  periodic  curve 


x(e)  =  e  +  b  sin  e 


y(e)  =  0.4  sin  e. 


(3.5) 


For  b  <1,  the  curve  is  a  single  valued  function  of  x. 

For  b  =  1,  the  curve  has  a  vertical  slope  at  e  =  it  , 
and  for  b>  1,  the  function  is  multivalued.  In  Fig.  3, 
the  curves  (3.5)  are  plotted  for  b  =  1,1.5,  and  2.0. 

The  map  function  u(e)  must  be  a  monotonically  increasing 
function  of  e.  Therefore  du/de  >  0  although  it  can  be 
exponentially  small  due  to  crowding.  The  functions  u(e) 
and  du/de  are  tabulated  for  the  curves  (3.5)  in 
Table  3.1.  Another  measure  of  the  crowding  is  given  by 
Jin  (g^)  .  In  Fig.  4,  £,n  (^)  is  plotted  for  various  values 

of  b  to  reveal  the  exponetial  nature  of  the  crowding 
phenomenon.  As  b  increases,  the  crowding  rapidly  becomes 
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severe  even  though  the  amplitude  of  the  wave  (3.5)  is 
quite  modest.  Similar  crowding  should  be  expected  in 
any  dynamic  simulation  of  a  breaking  wave. 


It  is  also  possible  to  formulate  a  set  of  differential 
equations  based  on  the  Menikof f-Zemach  approach  to  map 
a  time-dependent  boundary.  For  parametrized  boundaries 
of  the  form  (x (e, t) ,y (e , t) )  ,  the  mapping  function  u(e,t) 
is  determined  by 


3u 

3t 


3x  3_x 
at  3e 


(— ) 
'de' 


3y  3y 


at  3e 
2 


au 

ae 


(3.6) 


f  cot(!!.(e>-.u<e'l) 


-  If(e')  §*<«•) 


<af(e  >>2 


(|H(e'))2 


de ' 
2tt 


where  (g^)  ^  =  (-|^)  ^  +  (|^-)  ^  and  Ug(t)  is  chosen  so  that 
u(0)  =  0. 

Given  the  values  of  u,<j>,x,  and  y  at  some  time  t 
the  time  stepping  algorithm  proceeds  as  follows:  First, 
the  values  of  u(e),  x(e)  and  y(e)  are  used  to  determine 


the  map  derivative  du/de.  Next,  the  normal  velocities 
3<f>/3n  can  be  computed  from  (3.4).  Once  34>/9s  and 
3<J>/3n  are  known,  the  boundary  curve  (x  (e ,  t)  ,y  (e ,  t)  ) 
can  be  marched  to  the  next  time  step.  Then,  Bernoulli's 
equation  (1.2)  gives  the  boundary  values  of  4*  at  the 
next  step.  Finally,  the  map  Eq.  (3.6)  is  used  to  march 
u  forward  in  time. 

We  have  tested  the  time  dependent  mapping  equation  (3.6) 
on  the  mapping  of  the  region  bounded  by  a  cosine  curve  of 
increasing  amplitude. 


x(e,t)  =  e 

y (e,t)  =  t  cos (e) 


(3.7) 


and  on  the  regions  bounded  by  a  time  dependent  version  of 
the  breaking  curves  (3.5), 


x (e, t)  =  e  +  t  sin  e 


y  (e, t)  =  0.4  sin  e. 


(3.8) 


A 

A  fourth  order  Adams-Moulton  predictor-corrector  scheme 
was  used  to  march  the  map  function  u(e,t)  forward  in  time. 

At  the  times  tabulated  in  Tables  2  and  3,  the  mapping 
function  was  corrected  by  solving  (3.1).  The  time  integration 
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was  then  restarted  with  the  corrected  values  of  u(e). 

The  maximum  error  for  a  given  time  is  given  in  Tables 
2  and  3  for  32,  64  and  128  points.  The  minimum  of  the 
function  du/de  for  each  time  is  also  listed  to  give 
an  indication  of  the  crowding.  The  error  for  moderate 
distortions  was  fairly  insensitive  to  reductions  in  the 
time  step  At  but  was  reduced  markedly  when  the  number 
of  points  was  increased.  In  regions  of  severe  crowding 
the  time  step  must  be  very  small  in  order  to  ensure 
accuracy  for  an  explicit  integration  scheme.  Too  large 
a  time  step  can  destroy  the  raonotonicity  of  u(e). 

We  have  also  applied  the  integral  equation  (3.1) 
and  time-dependent  evolution  equation  (3.6)  to  the 
numerical  simulation  of  Rayleigh-Taylor  instability.  The 
initial  conditions  for  the  Rayleigh-Taylor  problem  are 
as  follows.  Fluid  of  density  1  lies  above  the  periodic 
interface 

y(e,t  =  0)  =  0.5  cos(e) 

x  (e ,  t  =  0)  =  e 

and  is  initially  at  rest.  Below  the  interface,  there  is 
a  vacuum.  The  resulting  flow  is  unstable  under  gravitational 
acceleration.  The  results  plotted  in  Fig.  5  are  obtained 
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using  the  integral  equation  (3.1).  With  60  points  per 
wavelength,  we  were  unable  to  continue  the  calculation 
past  a  time  of  t^3.5  at  which  the  amplitude  to  wavelength 
ratio  of  the  spike  (at  x  =  ±ir  )  is  about  5.4/2ir  -  0.86. 

The  degree  to  which  the  total  energy  and  the  rate  of  mass 
flux  are  conservedgives  a  good  indication  of  the  reliability 
of  the  simulation.  After  a  time  of  3.0,  there  is  a 
progressive  degradation  of  conservation  of  these  quantities. 
This  deterioration  is  also  reflected  in  the  spike 
acceleration.  For  large  t,  the  spike  should  be  nearly 
infreefall  with  an  acceleration  of  -1.0  in  our  units.1 
In  contrast,  the  present  simulation  shows  a  spike  acceleration 
which  decreases  (in  absolute  value)  below  1.0  after  t  =  3.0. 
Hence  we  conclude  that  the  results  are  not  reliable  beyond 
t  =  3.0.  Similarly,  the  time-dependent  evolution  equation 
(3.6)  gives  results  for  this  problem  that  are  reliable 
only  until  t;  3.0. 

The  present  conformal  mapping  methods  give  results 
for  Rayleigh-Taylor  instability  that  are  quite  good.  The 
amplitude/wavelength  ratio  has  increased  by  about  a  factor 
10  before  the  60  point  calculations  degrade.  Menikoff 
and  Zemach  (private  communication)  obtain  similar 
amplifications  before  their  calculations  break  down.  However, 
the  reasons  for  degradation  of  the  calculations  at  large  time 
remain  unclear.  On  the  one  hand,  the  conformal  mapping 
methods  described  in  this  Section  are  capable  of  resolving 
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1 


much  more  highly  deformed  interfaces  than  achieved  at 

breakdown,  even  with  60  points.  On  the  other  hand, 

1  2 

new  vortex  methods  '  have  been  used  to  calculate 
Rayleigh-Taylor  instability  with  similar  spatial 
resolution  to  at  least  twice  the  amplifications  achieved 
here.  It  seems  that  our  method  of  coupling  free-surface 
dynamics  and  conformal  mapping  introduces  numerical 
inaccuracies  (observed  as  rapid  oscillations  of  <j)  and 
0  for  t>  3.0).  It  is  possible  that  this  deficiency 
may  be  corrected  by  more  sophisticated  conformal  mapping 

3 

techniques  . 
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Table  1.  Mapping  functions  for  the  'breaking'  curves  (3.5) 


b  =  1 

.  0 

b  =  1. 

_5 

b  = 

o 

• 

CM 

u(  e) 

du 

u  (e) 

du 

u  (e) 

du 

e 

de 

de 

de 

0 

o 

• 

o 

1.9239 

0.0 

2.3575 

© 

• 

o 

2.7670 

IT 

4 

1.5706 

1.9648 

1.8956 

2.3030 

2.2126 

2.6351 

TT 

O 

2.9331 

1.4326 

3.4076 

1.4615 

3.8677 

1.4772 

3it 

4 

3.7681 

7.0083X10"1 

4.1579 

4.9413X10-1 

4.5352 

3.0688xl0-1 

TT 

4.0981 

2.0146X10"1 

4.3291 

4.7511xl0~2 

4.6125 

4 . 5852xl0~3 

5tt 

4 

4.1841 

7.7401xl0_1 

4.3387 

6.4110X10-4 

4.6131 

9.8707xl0-7 

3tt 

2 

4.4406 

4.2459xl0_1 

4.3440 

4.0460xl0~2 

4.6131 

3.2512xl0~5 

7tt 

4 

5.2656 

1.2741 

4.7624 

1.3178 

4.6920 

7.0508xl0_1 

* 

The  results 

were  checked  for 

32  and 

64  points  and 

agreed 

to  the 

5  significant  digits  given  here. 


Table  3.  Error  in  the  conformal  mapping  of  the  time-dependent 

* 

breaking'  curves  (3.8)  . 


t 

Maximum  Error  (percent) 

N  =  32 

N  =  64 

N  =  128 

.8 

1. 85*10-1 

4.5xl0~5 

-10 

3.2x10 

1.6xl0~10 

1.0 

7.48xl0-2 

2.0xl0“3 

3.1xl0-7 

3.8xl0-10 

1.2 

1.63xl0“2 

3.9xl0“2 

9.0xl0_5 

6.4xl0-10 

1.4 

1.99xl0-3 

— 

4.9xl0_3 

8.5xl0-7 

1.6 

1.02xl0-3 

_  _  _ 

9.5xl0-5 

*  The  time  step  is  At  =  0.001  and  48-bit  mantissa  arithmetic  is  used. 
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Figure  Captions 


Figure  1. 


Figure  2. 


Figure  3. 


Figure  4. 


A  schematic  plot  indicating  the  sequence  of 
conformal  maps  used  to  solve  inviscid  free 
surface  flow  problems.  Here  the  fluid  lies 
below  the  interface  y  =  n(x,t)  as  in  the 
water  wave  problem. 

A  plot  of  the  Stokes  wave  profile  at  t  =  0 
and  at  t  =  2t r.  The  amplitude  is  80%  of  the 
maximum  Stokes  wave  amplitude.  The  FFT 
time-dependent  mapping  equation  (2.21) 
is  used  with  N  =  64  points.  The  dots 
indicate  the  numerically  computed  position 
of  the  interface.  The  solid  line  is  obtained 
from  Pade  summation  of  the  perturbation  series 
for  Stokes  waves. 

A  plot  of  y  vs  x  for  the  'breaking'  curves 
(3.5),  x  =  e  +  b  sine,  y  =  .4  sine,  for 
b  =  1.0,  1.5,  2.0. 

A  plot  of  Jin  du/de  for  the  breaking  curves 
plotted  in  Fig.  3  (a)  b  —  1.0,  (b)  b  =  1.5 
(c)  b  =  2.0.  Here  the  Menikof f-Zemach  equation 
(3.1)  is  solved  for  the  conformal  mapping  function 
u(e).  Observe  the  exponentially  strong  crowding 
for  b  >  1. 
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Figure  5.  A  plot  of  the  interface  y(x,t)  for  the  Rayleigh- 
Taylor  instability  with  initial  conditions 
y(x,t  =  0)  =  0.5  cosx<Mx,t)  =  0  for  t  =  0.5 
to  t  =  3.5  in  steps  of  0.5.  Here  60  points 
per  wavelength  are  used .  Both  the  integral 
equation  (3.1)  and  time-dependent  equation 
(3.6)  degrade  significantly  in  accuracy  for 
t  >  3  at  this  spatial  resolution. 
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